Chaos in the Kepler System 



C. Chicone* , B. Mashhoont and D. G. Retzloff* 

The long-term dynamical evolution of a Keplerian binary orbit due to the emission and absorption 
of gravitational radiation is investigated. This work extends our previous results on transient chaos 
in the planar case to the three dimensional Kepler system. Specifically, we consider the nonlinear 
evolution of the relative orbit due to gravitational radiation damping as well as external gravitational 
radiation that is obliquely incident on the initial orbital plane. The variation of orbital inclination, 
especially during resonance capture, turns out to be very sensitive to the initial conditions. Moreover, 
we discuss the novel phenomenon of chaotic transition. 
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Q\ , I. INTRODUCTION 

0^ ■ 

' In a recent investigation of the radiative perturbations of the planar Kepler problem Q, we have found evidence 

Ph for transient chaos. Specifically, we have considered in our previous work an isolated binary system consisting of 
^ ^ , masses mi and m2 whose internal structures have been ignored, i.e. they have been treated essentially as Newtonian 
point masses with positions Xi and X2 in an inertial frame. The relative motion is damped since orbital energy and 
, angular momentum leave the system due to the emission of gravitational radiation. Moreover, the binary system 
■ is perturbed by an external normally incident plane gravitational wave. All radiative perturbations of the binary 
I I system have been treated in the quadrupole approximation in our work [Q; therefore, it follows that the center-of- 
^ r. mass motion is unaffected by the radiative effects. The well-known connection between the Kepler system and the 
f — ' harmonic oscillator implies that the system under consideration is analogous to a certain damped oscillator with 
. external periodic forcing. In fact, it is expected that replacing the emission and absorption of gravitational radiation 
by other similar damping and forcing mechanisms, respectively, might lead to similar phenomena. For instance, the 
influence of the external radiation could be replaced by the tidal perturbations of a distant third mass. 

The transient chaos in the dynamical evolution of the relative orbit under radiative perturbations appears to be 
' associated with the phenomenon of capture into resonance as explained in detail in our present work. Moreover, 
I we point out in this paper an interesting aspect of transient chaos, namely, the phenomenon of chaotic transition. 
The slow evolution of the binary orbit while locked in resonance can be explained using the second order averaged 
I , dynamics that is described here for the three dimensional Kepler system. 

^Jpj' It is important to emphasize the idealized nature of the model under investigation here We consider the 

. . simplest Keplerian model in which the effects of emission and absorption of gravitational radiation are taken into 
account. The radiative damping is of particular interest in view of the timing observations of the Hulse- Taylor 
relativistic binary pulsar PSR B1913-I-16: the data can be interpreted in terms of energy loss due to the emission of 
gravitational radiation in agreement with general relativity ||^,^ . For recent discussions of the measurable relativistic 
d ' effects in binary systems, we refer to the investigations of Kopcikin and his collaborators [Q. 

It is estimated that about half of all stars are members of binary or multiple systems; therefore, it is possible 
that the general phenomena we have encountered in our theoretical investigations have actually occurred in nature. 
In particular, our results should be relevant for the behavior of relativistic binary pulsars; indeed, a number 

of such systems have been discovered since the existence of the original binary pulsar PSR B1913-I-16 was first 
recognized [|[D. This possibility provides the impetus to analyze the Kepler system under more general conditions, 
especially in connection with the nature of its chaotic behavior. To this end, we consider in this paper the radiative 
perturbations of the three dimensional Kepler system. As work continues on relativistic binary systems — following the 
Taylor-Hulse work on the original binary pulsar — and more are discovered and studied, it may be that the theoretical 
results presented in our work could be useful in the elucidation of observed phenomena in such interesting astronomical 
systems. 
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II. DYNAMICAL EQUATIONS 



The simplest equation of relative motion for the Kepler system including radiative perturbations is given by 

'^"'^ -k^^+K =-eJC,,{t)x\ (1) 



where x :— xi — X2 is the vector describing relative motion, k :— GqM , M :— mi +m2, the quantity —TZ is the relative 
acceleration caused by gravitational radiation reaction, and e/C is the tidal matrix of the external wave evaluated at the 
center of mass of the binary system that we take to be the origin of spatial coordinates in the background inertial frame. 
Thus, in this frame the motion of mi and m2 can be described by xi(t) ~ (m2/M)x(t) and X2(t) = — (mi/Af)x(t), 
respectively. Various properties of the system have been discussed in detail in our recent investigations of the 
planar Kepler problem therefore, we concentrate in this paper on the long term deviations of the system (|^) 

from planar motion. 

In the quadrupole approximation, the standard expression for radiation damping may be reduced — for the Kepler 
system under consideration — to 

5°5|,|3 [V-^-x(x-v)x] , (2) 
where v — dx/dt is the relative velocity, x = x/|x|, and 

Ah Ah 

iP = Uv^ - 30(x . v)2 - — , x = 36i;2 - 50(i . v)^ + — . (3) 
|x| 3|x| 

Moreover, the tidal matrix for a plane monochromatic wave of frequency fl propagating in the (x^, a;^) -plane along a 
direction that makes an angle 0, < < tt, with respect to the z'^-axis is given by (cf. appendix A) 

/Cii = afl^ cos^ Q cos fit, 
/Ci2 = l3il'^ cos e cos{nt + 0o), 
/Ci3 = —afl^ cos 9 sin Q cos ilt, 
K,22 = — ail^ cos fii, 

/C23 = -/31^^sinecos(m + 0o), (4) 

and the other components of the tidal matrix are determined by the fact that (/Cy) is symmetric and traceless. Here a 
and /3 are amplitudes of the two independent linear polarization states of the incident radiation and 00 is the constant 
phase difference between them. The background inertial frame is fixed by the form of the tidal matrix and the fact 
that the initial unperturbed orbit is assumed to be in the (x^, a;^)-plane. The strength of the external perturbation 
is given by e, < e << 1, so that a and (i are of the order of unity; at present, it is expected that e ^ 10~^° though 
gravitational waves have not yet been detected in the laboratory. 

It must be emphasized that the simple form of the external perturbation assumed here is meant to represent the 
dominant component of a wave packet composed mostly of wavelengths much larger than the size of the Keplerian 
system. In the absence of the external perturbation, an orbit in the (x^, a;^)-plane would remain confined to this plane 
under radiative damping. If the external wave is turned on at a certain instant, it is expected that the equation of 
motion would represent the steady-state situation after transients have died away. 

It is interesting to subject the equation of motion to certain scale transformations that would render it dimensionless. 
To this end, let x = _Ro x and t — Tgi such that Rq and Tq are constants and a tilde denotes a dimensionless quantity. 
In this way the equation of motion keeps its form except that k must be replaced by k, 

k = > (5) 

U, = ft/To and the strength of the gravitational radiation reaction is given by the dimensionless quantity S, 

^ _ 4G'omim2 
" bc^ToRo ■ 



(6) 



In the rest of this paper, we fix the relationship between Rq and Tq by /cTq — Rq, so that fc = 1. Thus if Rq is taken 
to be the semimajor axis of the initial Keplerian orbit, then its period is given by 2ttTq. We adopt this convention 
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regarding the initial orbit throughout this work. It follows that S < (20-^2)"^ in all physically significant cases; in 
fact, S ^ 10^^^ for the binary pulsar PSR B1913 + 16 Henceforth we deal with dimcnsionlcss equations of 

motion in which all tildes are dropped. 

It is appropriate to introduce cylindrical coordinates {p,9,z) such that — pcos.9, = psin0, and = z. The 
equations of motion can then be expressed in the form 
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where we use equation (^) for the external perturbation and 

/C33 = -(/Cii +/C22) = sin^ecosllt. 
In equation (Q), the radiation reaction terms are given by 
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where ip and x ^re given by equation 

V- - 12 (pI + 



(p2 + ^2)3/2 

Sijjpe 

(p2 + 22)3/2 

__5__ 

(p2 + ^2)3/2 



i'Pz - X 



z{ppp + zpz) 



p-- 

and can be expressed as 

{PPp + zpzf 



P. 



30- 



X^i&\pl + ^+P. 



50 



p2 + z2 
p2 + z2 



(p2 



z2)l/2^ 



3(p2 +^2)1/2- 



It now remains to integrate these equations with appropriate boundary conditions. 
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III. NUMERICAL EXPERIMENTS 



The equations of motion (|^) have been so formulated that for 8 = 0, we recover the results of our recent work 
We wish to investigate how our previous results vary with Q as the direction of incidence of the external radiation 
deviates from the normal. It is clear that our previous theoretical results, such as regarding the persistence of periodic 
orbits, would still hold for sufficiently small O, < 1; however, the analysis for arbitrary O would be cumbersome 
and we therefore resort to numerical experiments. To express our numerical results in a form that could be easily 
compared with previous results |^,^, we note that the energy E, total angular momentum G and the z-component of 
the angular momentum H for the osculating ellipse in our system can be expressed in terms of cylindrical coordinates 
as 
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(11) 

(12) 
(13) 
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respectively. The corresponding Delaunay elements are the action variables {L,G,H), where L = (— 2£')~^/^. The 
osculating ellipse is depicted in figure |^. 

Some typical results of our numerical experiments are presented in figures |[^- Figure |^ explores the phenomenon 
associated with capture into resonance; here the osculating orbit has eccentricity e, G = L(\ — e^)^/^, and H — Gcosi, 
where i denotes the orbital inclination. The angle of incidence is = tt/IO. Figures |^-|^ explore the phenomena 
associated with transient chaos for O = 7r/6. 

To arrive at these representative results, we begin our numerical work with a resonant orbit given at some initial 
instant of time by 

{pp,pe,pz,p,0,z) = (eo, 1,0, 1,0,0), 

which corresponds to an osculating ellipse in the (x^, a;^)-plane with eccentricity gq, unit orbital angular momentum 
about the z-axis, semimajor axis (1 — Gq)^^ and Keplerian frequency (1 — Bq)'^/^. In figure ^ for instance, this initial 
orbit would have h = 0, g = ^-k /2 and ij — 7r/2. We then assume an external wave that is a coherent superposition 
of the two independent linear polarization states with equal amplitudes (a = (3) and zero phase shift between them 
((/iQ — 0). The wave is incident at an angle 8 that we vary in our experiments. We choose the wave frequency 17 = 
to correspond to the (m : 1) resonance. Integrating forward and backward from this initial system, we explore the 
region around the resonance configuration. Starting from a state before resonance capture, we then integrate forward 
until the integration routine gives up due to the rapid collapse of the binary system. To compare our results with 
the planar case = 0, we choose the two planar configurations with initial Cq — 1/2 that we studied in a previous 
paper ||l|. In the first system, i.e. figure 1 of ||l[], we used e — 10"**, 5 — 10"^, a — 13 = 2, = Q and O = L^^ with 
io = (1 — Gq)^^/^. This (1 : 1) resonance is studied for 6 = tt/IO in figure || of the present paper. In the second 
system, i.e. figure 2 of H, we used e = lO'^, 5 = 10-^ a = /3 = 2, 0o = and = 2L-^^ with Lq = (l-eg)-!/^. This 
(2 : 1) resonance is studied for 8 = 7r/6 in figures of the present paper; in these figures, every tenth iterate of the 
2-11/0, stroboscopic Poincare map is plotted. This second case exhibits transient chaos and to be reasonably certain 
of the various features of this chaotic regime, we have used two different standard routines for stiff integration of (Q) 
using double precision arithmetic. The results of the integration of this configuration are depicted in figures ^ and^. 
These are to be compared with figures ^ and ^ respectively, which depict the results of integration of essentially the 
same system albeit with a different starting point. It is important to note that the orbital inclination is particularly 
sensitive to transient chaos as is clear from the behavior of cos i in figures ^ and ^. Moreover, as the binary system 
collapses, we expect that eventually the binary orbit would tend to a circle (e 0). In figures the endpoint 
of integration exhibits a sharp drop in eccentricity as would be theoretically expected. Note that this behaviour of 
eccentricity coincides with the sharp bend in L{t) near L = 1 corresponding to chaotic transition as described in 
section 

In the following sections, we present a theoretical explanation of some of the main features of the behavior of the 
system over a long period of time. For this purpose, essential use is made of the dynamical equations in terms of 
Delaunay elements (cf. appendix B). 



IV. AVERAGING 



The long-term behavior of the dynamical system under investigation is best revealed after we average the system 
over the Keplerian frequency lo = L~^. In fact, the Delaunay equations (cf. appendix B) show clearly that the three 
dimensional Kepler system has only one intrinsic frequency, namely, lo. Let us write the equations of motion in the 
form 
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(14) 



where Ti* is the Hamihonian associated with the external perturbation, A = 5/e and /_d, D e {L,G,H,£,g,h}, 
denotes the gravitational radiation damping. Here t — ilt is a new angular variable. 

In the absence of the external perturbation, we have shown Q that the system maintains its orbital plane — i.e., the 
inclination of the orbit remains fixed — but loses energy and angular momentum so that the orbit decays monotonically 
and eventually collapses. During this process the orbit tends to a circle, i.e. the eccentricity e decreases. With the 
external radiation present, as in (|^), we have a dynamical system with two frequencies ui and fl and the possibility 
of resonance needs to be taken into account. This would occur if relatively prime integers m and n exist such that 
muj = nfl; the (to : n) resonance manifold is then given by the set 



{iL,G,HJ,g,h,T) 



m— = nfl}. 



Let us first assume that the dynamical system (^_4|) is off resonance; then, averaging (|14 | ) o ver a period of the external 
radiation results in a system that is simply damped. That is, the terms involving Ti.* iii([l4|) simply drop out; therefore, 
the system tends to a circular orbit as it decays on the average while maintaining the inclination of the orbital plane 
(cf. appendix B). This explains figure |^ before and after resonance capture; that is, in these regimes e decreases 
monotonically on the average while cosz is essentially constant. If the resonance condition is satisfied, then the 
system either passes through the resonance without being captured or is captured into the resonance. In the former 
case, the system behaves on the average in the same way that it would off resonance. In the latter case, however, the 
resonance condition fixes the semimajor axis of the osculating ellipse on the average. Thus the external wave deposits 
energy into the system so as to balance radiation damping on the average and hence L remains fixed on the average 
while the exchange of angular momentum takes place. The resulting change in G and H affects the shape as well 
as the configuration of the orbit; in fact, the eccentricity and the inclination of the orbit vary while its semimajor 
axis oscillates about the value fixed by the resonance. This is clearly indicated in figure H, where the inclination and 
the eccentricity of the orbit both decrease as a result of passage through resonance. The rate of damping depends 
on the shape of the orbit and this fact accounts for the difference in the slope of L{t) before and after capture into 
resonance. During resonance capture, we let L — Lq + t^^^T) and I = (f> + nflt/m, where Lq^ — nVt/m and T> and (j) 
are the new canonical variables associated with resonance. The behavior of the dynamical system while trapped in 
resonance can be studied using the method of partial averaging (cf. appendix C). It turns out that for the radiative 
perturbations of the three dimensional Kepler system under consideration here resonance is possible only for n = 1. 
Numerical experiments then reveal that the (1:1) resonance has a simple structure as in figure 2. More complicated 
structures are expected for higher order resonances as described in our previous planar work In the present paper, 
we explore a (2 : 1) resonance in figures 3-6. A general discussion of the (to : 1) resonance for the Kepler system is 
beyond the scope of this work. 

To explain the behavior of the system while trapped in resonance in a quantitative manner, we need to develop 
the first and second order partially averaged dynamics for (^^ following the method adapted in our recent work on 
the planar Kepler system |2] . This is done in appendix C and the resulting system of second order partially averaged 
Delaunay equations for the (to : 1) resonance are 
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where 7^ and Tg are coniplicated expressions involving a, /9, fi, (/)o, O and the Delaunay elements given by (C25) in 
appendix C. It is important to note that the first-order averaged equations are formally the same as in |^]; that is, to 
order e^/^ the orbital inclination is fixed and we simply have an oscillator in V and a pendulum with constant torque 
in (j) as before. To second order in e^^^, the damping (or antidamping) term enters the pendulum equation and this 
results in an antidamped (or damped) harmonic oscillator in T) as in the case illustrated in figure 2. These features 
have been explained in detail in therefore, it is more interesting to investigate the new aspects of resonance capture 
associated with the slow oscillation of the plane of the osculating orbit. In fact, the possibility of the variation of the 
inclination to second order in e^^^ would result in an additional feature that is required for a complete explanation of 
figure]^. The pendulum in (p will couple, in the general case, to an equation for the variation in the orbital inclination. 
To develop this latter equation, we note that H = Gcosi and hence d{cosi)/dt = {HG — HG)/G^, where G and H 
are given by (p^. We are interested in the slow variation in the inclination of the orbit during resonance capture. The 
average behavior of the orbital inclination is obtained if we substitute for G and H from the second order partially 
averaged system ( p^ ) instead. In this way, we find 
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(16) 



for the slow variation of the orbital inclination while the system is captured in resonance. It is remarkable that the 
explicit damping terms proportional to A in ( p^ ) drop out in equation (|l6|); indeed, in the absence of external waves 
there would be no resonance capture and cosi would be constant. Thus the rate of variation of cosi is directly 
proportional to the presence of the external perturbation as in ([l6| ) . It is difficult to draw exact conclusions regarding 
cosi from ( p^ ) and (p^); however, we intuitively expect that the orbital inclination would undergo on average simple 
oscillatory movements during resonance capture. These oscillations are expected to be slow, since cos i varies in time 
over the temporal scale given by et. These conclusions are generally consistent with the results of our numerical 
experiments. 



V. TRANSIENT CHAOS 



Numerical experiments suggest that our system exhibits chaotic transients near the "exit from resonance" . In 
our experiments, we considered the (2:1) resonance similar to the planar case in figure 2 of ||l|]. Starting at the 
point {pp,pe,Pz, P, Q, z) = (0.5, 1, 0, 1, 0, 0) which is near the resonance manifold, backward integration shows that the 
system is captured into resonance; then, after a sufficiently long forward integration, the orbit leaves the vicinity of the 
resonance and continues toward collapse. However, the length of time required for the system to exit the resonance is 
very sensitive to the integration method and to changes in the initial conditions. Also, plots of L versus t show that 
the signal appears to pass through a chaotic region after it finally exits the resonance. This strongly suggests that 
the system is passing through an unstable chaotic set; that is, there is transient chaos. 

We expect transient chaotic effects to be present throughout the resonance capture region; however, in our ex- 
periments transient chaos appears most prominently after exit from resonance. Following the exit from the (2:1) 
resonance under consideration, the system collapses relatively slowly until it undergoes a certain chaotic transition 
and the relative orbit that emerges after this transition collapses much more rapidly. In fact, the system makes a 
transition from one relative orbit to a totally different relative orbit in a rather short period of time while undergoing 
what appears to be transient chaos. For lack of a better characterization of this unexpected phenomenon, we refer 
to it as "chaotic transition" . As a consequence of this transformation, a "bent knee" appears in figures near 
L = 1; indeed, the rate at which the system collapses suddenly changes at the "bend". That is, the orbital eccen- 
tricity increases very rapidly and this leads to a more rapid rate of loss of energy due to the emission of gravitational 
radiation — since this is proportional to (1 — e^)^*"/^ as can be seen from inspection of equation (15) of ||l| — and hence 
the "bend." We have verified by two different numerical methods that these chaotic phenomena are stable features 
of the (2 : 1) resonance. 



VI. CONCLUSION 



We consider in this paper a three dimensional Kepler system that is subject to radiative perturbations caused by 
radiation damping as well as an incident monochromatic gravitational wave. 

A Keplerian binary system constantly loses energy to gravitational radiation according to general relativity. There is 
a lack of complete reciprocity between the emission and absorption of gravitational radiation, however. The absorption 
of gravitational radiation by the binary is not monotonic and the system sometimes gains energy from the external 
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wave and sometimes deposits energy into the wave. The behavior of the system averaged over the period of the 
external wave is thus one of continuous coUapse due to radiation damping except when it is captured into resonance. 
Other than at resonance, the average behavior of the system is that the orbital plane remains fixed while the system 
shrinks as it tends to a circle (a — > 0,e — > 0). At resonance, the system on average steadily gains energy from 
the external wave in order to balance the steady loss of energy to radiation. The resonance condition muj — nfl is 
necessary but not sufficient for the occurrence of this delicate balance at the (m : n) resonance. In general, transient 
chaotic behavior is expected near a resonance. As the system evolves, it passes through and is affected by a dense set 
of (m : n) resonances. When it is indeed captured into a resonance, we find that n = 1. 

We have not proven the existence of transient chaos in the Kepler system when radiative perturbations are taken 
into account. Instead, we have presented further numerical evidence in support of this hypothesis in this paper. Two 
different codes for the numerical integration of equation have been employed and they have yielded qualitatively 
the same results. Our conclusions regarding transient chaos and chaotic transition are based on these numerical 
results for the three dimensional Kepler system. 

The parameter space of the problem under consideration is very large. No significant attempt has been made at 
exploring this parameter space; for instance, we cannot exclude the possibility of existence of a strange attractor 
in our system. Instead, we have concentrated our attention on a (1 : 1) resonance and a (2 : 1) resonance in the 
three dimensional case that are familiar from our previous planar work The (1:1) resonance has a simple 

structure and its general features can be explained by the second order averaged dynamics. That is, the semimajor 
axis undergoes antidamped oscillations about its resonance value while the orbital inclination slowly oscillates until 
the system leaves the resonance; on the other hand, angular momentum is transferred to the binary orbit during 
resonance such that the orbit after the resonance is substantially different in eccentricity from the orbit before it 
gets captured into resonance and hence so is the rate of collapse of the binary due to the emission of gravitational 
radiation. 

The (2:1) resonance has a rich structure and much of our numerical analysis has been concerned with its elucidation. 
It appears that this structure is associated with transient chaos. Once the system leaves the resonance, the general 
trend towards collapse is accompanied by complex structure that we characterize as transient chaotic behavior. This 
is followed by a peculiar chaotic transition which may be another manifestation of transient chaos. The result is that 
the system collapses extremely rapidly as a result of going through this transition. Integrating backward in time from 
the resonance, the system grows in size. Since the radiation reaction force decreases with size at least as while 
the tidal force of the external wave grows as r, the end result is that the radiation reaction can become effectively 
negligible. This is the regime of Arnold diffusion in our problem and can lead to gravitational ionization as described 
in detail in our previous work |3[^. 

It is possible that for the general (m : 1) resonance other interesting phenomena may occur that are beyond the 
scope of this investigation. 

APPENDIX A: ABSORPTION OF GRAVITATIONAL WAVES 

Consider a plane monochromatic gravitational wave with wave vector K and frequency Q — c|K|. The spacetime 
metric is given by (/^^ = 7]^^, + e/i^,y, where we impose the gauge condition ho^ = 0. Here ?7^t^ is the Minkowski metric 
with signature +2; moreover, it is important to note that our gauge has the characteristic property that the worldline 
of a static test particle is a geodesic and this is consistent with the binary system as a whole remaining at rest in 
this external radiation field. The gravitational wave amplitude is then characterized by a symmetric and traceless hij 
with djhij = and Ohij = 0. Let 

hij = Re [hij exp{—int + iK. ■ x)], (Al) 

where hijK^ ~ 0. We are interested in the gravitational influence of this wave on a Keplerian binary system. If the 
wavelength of the radiation is much larger than the dimension of the system, then the interaction of the wave with 
the system is predominantly tidal in character. Then the tidal matrix is given by 

^^.=-^^(i,0), (A2) 

since the center of mass of the binary system is fixed at the origin of the spatial coordinates. We choose the coordinate 
system such that the wave is incident in the (a;, z) -plane; hence, K = (sin 0, 0, cosO) is the unit propagation vector 
of the wave. Let us then choose (K, N,y) to be an orthonormal triad with N — (cos 0, 0, — sinO). Expanding 
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the symmetric and traceless hij in terms of this basis triad and taking the transversahty of the wave into account 
(hijK^ = 0), we find 

/ly = 2a{N,Nj - y,yj) + 2/3exp(-i</)o)(Af,yj- + N^m), (A3) 

where 2a and 2/3 are the real ampUtudes of the two independent hnear polarization states and 0o is the constant phase 
difference between them. Here the zero of time is so chosen that a is real. It is now possible to recover equation (^) 
from 

The absorption of tidal gravitational energy from such a wave is not unidirectional in general; that is, a system can 
gain or lose energy as a result of its tidal interaction with an incident gravitational wave. On the other hand, we have 
found that during resonance capture the system always gains energy. 

The nonlinear evolution of the binary system has been emphasized in our recent work when the external perturbation 
can, in effect, be considered to be a plane monochromatic gravitational wave However, we have not considered 

an arbitrary incident wave packet or a stochastic background of gravitational waves. These cases have been treated 
in the absence of radiation damping for the linear evolution of the binary system |^ . It turns out that for a binary 
immersed in a random, isotropic and unpolarized linear gravitational wave background characterized solely by its 
energy spectral density, the osculating elements of the relative orbit undergo random walks The extension of 
these results to the long term nonlinear evolution of the system remains a task for the future. In this connection, let 
us mention that nonlinear stochastic differential equations have been the subject of numerous investigations (see, for 
example, 



APPENDIX B: EQUATIONS OF MOTION IN DELAUNAY'S VARIABLES 



The equations of motion in terms of Delaunay's action-angle variables are of basic importance for the dynamical 
considerations in this paper. The standard treatments usually involve Hamiltonian systems; therefore, we present 
here a direct derivation of these equations since our system is dissipative. 

The equations of motion for the perturbed Kepler problem ([l|) can be written in vector form as 



\ = F 



(Bl) 



where r = |x| and F is given by 
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At any instant of time i, x and v specify the instantaneous orbital plane that contains the osculating ellipse with its 
focus always at the origin of coordinates. Let us choose the coordinate system such that the osculating ellipse at time 
t has the form given in figure ^ 

To obtain the equations of motion in Delaunay elements, we employ an orthonormal frame field adapted to the 
instantaneous osculating ellipse, write the component equations of motion using this basis, and perform a transfor- 
mation from the inertial Cartesian coordinates to those adapted to the instantaneous osculating ellipse. We define 
the frame field as follows. Let f be a unit radial vector, x = rf , and n be a unit vector normal to the plane of the 
instantaneous osculating ellipse. To complete this frame field, we include s = n x f which is a unit vector in the 
orbital plane normal to f . Then F can be written as 



F = F^f + FsS + F„n 



(B3) 



using its radial, sideways and normal components. In what follows, we use the definition of the Delaunay elements 
(L, G, H, £, g, h) given by 

L = a^'\ 
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H = Gcos i, 
i = u — e sin u, 

g ~ argument of the periastron, 

h = longitude of the ascending node, 



(B4) 



where a is the semimajor axis of the osculating ellipse, e is its eccentricity, u is the eccentric anomaly, and £ is the 
mean anomaly. Only positive square roots are considered throughout this paper. The equation for the radial position 
r in terms of the true anomaly v is given by 



a(l - e2) 



The orbital energy IlLWj is given by 



where 



E = -v^ - - = , 

2 r 2a 



a(l-e2) 



(B5) 



(B6) 



using the fact that G = |x x v|. It follows from these relations and 

esinv 



It is clear from (Bl) that 



G 



dE 

— = F • V. 

dt 



It then follows from (B7), (B8), and the definition of L that we have 

Fre sin?) + 



— = F • V 

dt 



1) that 



a(l - 



(1 - e2)i/2 

To obtain the equations of motion for the remaining Delaunay elements, we use the relations 

d 



(B7) 



(B8) 



(B9) 



to get 



and 



Gfi, 



dG 
Itt 



dt 



(x X v) 



X X F 



G^ = s. 

dt 



(BIO) 



(Bll) 



Now consider a coordinate system [x' , y' , z') oriented to the instantaneous osculating ellipse such that the perturbed 
motion is given by x' = rcos {v + g), y' — rsin (v + g), and z' — 0. The transformation from the coordinate system 
(x,y,z) to {x',y',z') consists of a rotation by an angle h about the z-axis followed by a rotation by the inclination 
angle i about the line of the ascending node, which is the x'-axis. This transformation can be written as 



(B12) 



cos h 


— sin h cos i 


sin h sin i \ 


'y'] 


sin h 


cos h cos i 


— cos/isini 







sini 


cosi / 





It follows that the relative position is given by 
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X = r sin cos ip = r [cos h cos (v + g) — sin /i cos i sin (-0 + g)] , 
y = r sin sin = r [sin /i cos (v + g) + cos h cos i sin (v + g)] , 
z = r cos d = r sin z sin (w + g) . 

In the (a;', y', z') coordinate system, n = (0, 0, 1). Thus in the inertial frame, 

n= (sin /isini, — cos /isini, cos i), 

and 



(B13) 



s = n X r 

= (— cos h sin {v + g) — sin h cos i cos {v + g), 

cosi cos h cos {v + g) — sin h sin (-0 + 5), sin i cos (i) + g))- 



(B14) 



Substituting (B14) into (Bll) results in 

di rFn 

_^_cos(. + g), 

dh . . rFn . , ^ 
— sm I = — - sm (u + ,g j . 

From the definition of iJ in (B4) and equations ( B10| ) and (B15), it follows that 

dH , . 

= r[Fs cosi — Fn sm? cos (v + gj\. 



(B15) 



(B16) 



The dynamical equation for the mean anomaly is obtained by exactly the same method as used in (cf. appendix 
B in [|l|). The result for the three dimensional case is obtained by replacing Fg in [|l| by i^s; hence, 



d£ 
dt 



e{a] 



Y7^[Fr{~2e + cosv + ecos v) — Fs{2 + ecosv) siiiv]. 



(B17) 



To obtain the equation for the dynamics of the argument of the pcriastron, consider the relative velocity vector in 
the form 



Project onto the z-axis to obtain 



V = rr H s. 

r 



G 



z = r sin z sin {v + g) A sini cos {v + g). 



Differentiating the z-component of (B13) with respect to time and equating the result to (B18) gives 

d{v + g) G rFn cos i 



dt 



G sini 



sin [v + g). 



Now, by taking the time derivative of the logarithm of (B5), we obtain 



_ G _ G 
dt e 



-Fr cos V + Fg \ 1 



a(l-e2) 



smv 



(B18) 



(B19) 



(B20) 



Subtracting (P320| ) from (B19) gives the dynamical equation for the argument of periastron as 

dg rFn cosi 



dt 



G sinz 
G 



sin [v + g) 



'Fr cosv + Fs 1 



a(l-e2) 



smv 



(B21) 



To summarize, the dynamical equations in Delaunay elements are 
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-- = rF, 



dL 

n 

dG 

dt 
dH_ 

~dt 

de 

di 
dg 
dt 



(l-e2)i/2 



FrB sinv + Ft 



a(l 



r[Fs cos i — Fn sin i cos {v + g)] , 

-Yyj[Fr(— 2e + cosw + ecos v) — Fs{2 + ecosv) sin-O], 



e(a)iA 
rFn cos i 

G sin i 

(a(l-e2))i/2 



sin [v + g) 

Fr cosv + Fs I 1 



a(l-e2) 



smw 



dh rFn sin (w + g) 



G 



Sim 



(B22) 



It now remains to express Fr, Fg and Fn in terms of Delaunay elements. 

The equations of motion in Delaunay's variables can be expressed in the form given in To this end, let 

1 . 



be the Hamiltonian associated with the external perturbation and S = eA. Then (B2) can be written as 

F = -eVn* + eAf , 

where 



16 



f = — I 2Av^ - 20f^ + -^\r- — [ 12v 



G 



- I s. 

r 



(B23) 



(B24) 



(B25) 



Thus fr — f(x ~ '4')/i'^, fs = —Gxjj/r^ and /„ = 0. It is important to note that the functions ip and x can be 
expressed solely in terms of Delaunay variables i n the plan e of t he orbit {L, G, £, g) just as in our previous work 
It then follows from the inspection of equations (B24) and (B25) that when radiation damping acts alone, the orbital 
pla ne re mains fixed in space while the orbit tends to a circle as it shrinks (e — *■ and a — > 0). Using these results 
in (B22), one can show that equation (|lj) is recovered with appropriate definitions of the damping functions fu and 
the replacement of ilt by r. 



APPENDIX C: PARTIAL AVERAGING 



This appendix is devoted to the determination of the resonant behavior of the three dimensional Kepler system 
when averaged over the binary period. Averaging over this "fast" motion should reveal the "slow" oscillatory behavior 
of the system while in resonance and its exit out of the resonance. 

The perturbed Kepler system continuously emits gravitational waves and hence loses energy and angular momentum 
to gravitational radiation damping. Hence the semimajor axis of the osculating ellipse must steadily decrease except 
when the system is trapped in a resonance. At resonance, a delicate balance exists between the external gravitational 
perturbation and the radiation damping since the energy (and hence the semimajor axis) of the osculating ellipse 
would be constant on the average. That is, during resonance capture the energy deposited into the orbit by the 
external wave would equal on average the energy lost via the emission of gravitational radiation. In general, the 
system slowly drifts out of the resonance when this subtle balance is significantly upset. 

Let L = L{) when the resonance condition muo — nil is exactly satisfied; then, we assume 

L = Lo + e^/^V, (CI) 
e^^t + <p, (C2) 
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during resonance capture. As before we recast the equations of motion in Delaunay variables (|lj) in terms of the 
new variables (I?, 0). To this end, 
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1_ 3^^1/2 + 6^e + 0(e3/2) 



(C3) 



and so we obtain from upon Taylor expansion 



V 



G — —eF22, 



H 



12, 



r4 ' "V r5 



9 = £-^52, 
h = 6^62, 

where terms of order higher than e have been neglected. Here (-Fy) are defined by 



(C4) 



while 



-F22 

^32 

F42 
F52 
^62 



dn* 

dLdl ^ 'dL 



A/l, 
dfL 



dn* 
an* 

~dh 



dn* 
~dr 
an* 

'dG 

an* 

~dH 



A/g 



A/h, 



A/,, 



A/„ 



AA. 



(C5) 



(C6) 



Note that the first index in Fij refers to the equation in which it appears while the second index refers to its order 



in powers of the perturbation parameter e^/^. Moreover, all the terms on the right sides of ( |C5[ ) and ( |C6|) are to 
be evaluated at {LQ,G,H,(f) + nil.t/m,g,h) so that (Fij) become functions of {G,H,(p + nnt/m,g,h,t), since n* is 
explicitly dependent upon time. Let us further note that at resonance the appropriate perturbation parameter turns 
out to be e^/^, since only with this choice would the conjugate variables 2? and (j) both vary predominantly on the 
same slow temporal scale given by e^^^t. 

The system ( |C4| ) is periodic in t with period 2nm/n. To average this system ov er i ts period, we introduce an 
averaging transformation that has the function of rendering the transformed system ( |C4D in a form that is averaged 
to first order in e^/^. The resulting system will then be replaced by the second order averaged system, where the 
terms of order e are simply averaged over the 27rm/f7 period. The averaging theorem ensures that the solution of the 
second order averaged system is sufhciently close to the solution of (C4) over a timescale of order e^^/^ ||^. 



Let us first define the average of Fij to be 



n 

27rm 



27rm/a 

F,j (G, H, — t + cj), g, h, t) dt. 



(C7) 



Since the averaging transformation involves only first order quantities, the only term of interest in ( |C4D would then 
be Fii and we let A(G, H, (f), g, h, t) = Fn — (Fu) be its deviation from the average. Next, we define A(G, H, (j), g, h, t) 
to be the antiderivative of A with respect to t such that (A) — 0. The averaging transformation is then given by 



V = V-e'^^AiG,H,c^,g,h,t), 



(C8) 
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G^G, H = H,(I) — (f), g = g and h = h. He nce the averaging transformation is on average equivalent to the identity 
transformation. Upon this transformation, (|C4|) takes the form 



v = 




G = 


-eF22, 


H = 




$ = 


_gl/2 ( A 






9 = 




h = 





3_dA 

Lq d(j) 



6 ^2 
^0 



3 A 



R 



42 



(C9) 



where terms of order e'^/^ and higher have been neglected. We note that (C9) is in averaged form to first order, but 
not to second order. Therefore, we now average (|C9| ) to obtain the second order partially averaged system that we 
wish to study. To simplify the notation, we write our main result in the form 

G^-e{F22), 
H = -e{F32), 

3 



9 = £(^^52) 

h = e{F62) 



V 



P2 + (F4; 



(CIO) 



since (A) = and {dk/dcj)) = 0. We remark that ( CIO ) is simply the averaged form of our original system (|CJ); 
however, this will not be true in general but happens to be the case here. It now remains to compute the averages 
(Fy ) defined by (|C7|). 

Let us first focus attention on the radiation reaction terms. It follows from the discussion in the last paragraph of 
appendix B that the effect of this frictional force is to maintain the orbital plane while the relative orbit loses energy 



and shrinks (a - 
and is given by 



0,e ^ 0). Therefore, the result of averaging is essentially the same as in our previous work 



(/l) 
(/g> 



1 



73 



37 



+ T^e + — e 



{fa) cosi. 



3 



12 



and (/,) - if,) = (/„> 



0. Moreover, {dfL/dL) ~ d{fL)/dL and hence 

1 



' dL' 



(146 + 37e^). 



(Cll) 

(C12) 
(C13) 



(C14) 



All these quantities are to be evaluated at L = Lq. 

We consider next the terms involving the external Hamiltonian eTi*, which must first be expressed in terms of the 
Delaunay elements. To this end, let 

S{h) — sin h, 
U{G, H, h) — cos cos i S + sin 6 sin i, 
V{G, H, h) — cos h cos i, 

M^(/i) = cose cos /i, (C15) 



and define Pa and Pa, cr — 0, ±, such that 
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P- = -SV-UW, 
Po = SW - UV, 
P+ = SW + UV, 

P- = -SU + VW. (C16) 



Using these quantities, we can form Va{G, H, h,t), 



= ]^an^ cosVlt P„ + ]^l3n^ cos{Vlt + (t)Q)P„. (C17) 
It is now possible to express Ti* in terms of the Delaunay variables as follows: 

(7 = 0, ± 

where 

Q°(L, G, £, <?) = = a^{\ + ^e^) - ^ ^ J,.(i^e) cos vi, (C19) 

1^=1 

5 °° 

G,£,g) = cos(2t; + 2^) = -a^e^ cos25r + ^(A,, cos2gcosi'^ - B,, sin 2^ sin z^^) , 

5 °° 

{L, G, e, g) = sin(2w + 2g) = -a^e^ sin 2^ + ^{A^ sin 2g cos i'^ + cos 2g sin vtj. 



(C20) 



Here A^, and are functions of the eccentricity e = [L? — G'^Y^'^ /L and can be expressed in terms of the Bessel 
function Jv(x) as 

A, = -^[2ve{l ~ e^)Jl(ve) - (2 - e^)J,{ve)], 
V e 

= -^(1 - e'f'^ieJliue) - v{l - e')Uue)], (C21) 

where J'y(x) = dJ^ix) / dx. 

We are interested in the average value of 7i*; therefore, it is convenient to express Vcr as 

Va = Cff cos l^i + Sa sin (C22) 



where C„ and 5o- can be explicitly determined using (C17). Moreover, it proves convenient to express the Fourier 
series expansions of Q'" collectively as 

oo 

= a'^ + Y^{aZcosvl^bZsmvl). (C23) 

The average of 7Y*, 

{n*)^- / n*{L,G,H,—t + <l>,g,h,t)dt, (C24) 

is given by 

{n*) = %{L, G, H, g, h) cos m<j> + % {L, G, H, g, h) smmcj), (C25) 
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where 



2 ^ 



(C26) 



for n = 1. If 71 7^ 1, then (7i*) = 0; hence, we defi ne 7^ = 7^ = m this case. 

It is now possible to express the quantities in (CIC) that depend on the external perturbation in terms of T^c and 
7^. For instance, it is straightforward to use the formulas (C18), (C23) and (C25) to prove that 



d_ 



{n* 



Moreover, the other partial derivatives of TC* with respect to the other Delaunay variables simply commute with the 
operation of averaging. For instance, {&H* /dL) = d{Ti*) / dL, which is finally evaluated at L — Lq. The second order 
partially averaged equations would thus have a form similar to that for the planar case except for additional equations 
for the slow variation of the extra Delaunay variables H and h. More explicitly, the partially averaged equations for 
the three dimensional Kepler system are given by equations (19) of (with Lq, </?—!■(/), Tc —> and Tg — > 7^) 

together with 



H = -e 
h 



dh 



cosmfi 



dh 



A 



sm mfl 



+ 7e^) cosi 



The full set of partially averaged Delaunay equations is given by (|15| 
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FIG. 1. Schematic plot of the osculating ellipse in the three dimensional Kepler problem. The instantaneous position vector 
is x(t) = (x, y, z). The unit vector n is normal to the instantaneous orbital plane and points in the same direction as the orbital 
angular momentum. 

FIG. 2. The plots are for system (^ with parameter values e = 10"'', S/e = 10~^, a = 2, l3 = 2, (po = and 6 = n/10. 
Here, Lo = (4/3)^''^ 1.1547 corresponds to (1 : 1) resonance with fl = 1/Lq. The top panel shows L = a^^^ versus time t for 
the initial conditions {pp,pg,p^, p,e,z) equal to (0.17554040552969, 0.57754077894430, -0.050059905721582, 2.7059982365422, 
—6.7506347480153, 0.14058165333520). The middle panel shows cosi versus time and the bottom panel shows e versus t. 



FIG. 3. The plots are for system (|7|) with the initial values of {pp,pg,pz, p,0, z) equal to 
(-0.6579651108853, 0.74664767379829, -0.12269473267506, 1.4254331835434, -10.1974793053511, -0.052493593111567) and 
parameter values e = 10"^, S/e = lO"'', a = 2, /? = 2, 00 = 0, 6 = n/6, and = 2/Lo with Lo = (4/3)^^^ so that the initial 
value is near the (2 : 1) resonance. The top panel depicts L versus time, the middle panel depicts G versus time, and the 
bottom panel depicts the eccentricity e versus time. Every 10th iterate of the 2-r/Q, stroboscopic Poincare map is plotted. 

FIG. 4. The top panel depicts cos i versus time, the middle panel depicts G versus L and the bottom panel depicts H versus 
L with the same data as in figure ^ 

FIG. 5. The plots are for system (Q) with parameters as in figure ^ but with initial values of {pp,pg,pz, p,6, z) equal to 
(-0.013456419026, 0.510820155326, 0.161909540107, 2.35369782619, -644207.714939, -2.17211570796). The top panel depicts 
L versus time, the middle panel depicts G versus time, and the bottom panel depicts the eccentricity e versus time. Every 10th 
iterate of the 2-k/^I stroboscopic Poincare map is plotted. 

FIG. 6. The top panel depicts cos i versus time, the middle panel depicts G versus L and the bottom panel depicts H versus 
L with the same data as in figure pi 
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